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Abstract —  A  spectral  domain  method  is  presented  to  calculate  the  potential  of  a  vertical 
dipole  in  a  multilayered  medium  as  a  model  of  long-range  radio  propagation.  The  spectral  domain 
Green’s  function  (SDGF)  for  structures  with  up  to  hundreds  of  thousands  of  layers  is  calculated 
using  an  efficient  matrix  formulation,  thus  enabling  simulation  of  continuously  stratified  media. 

The  SDGF  is  sampled  to  perform  pole/residue-extraction  using  contour  quadratures  and  to 
perform  a  novel  asymptotic  Filon-Clenshaw-Curtis  (FCC)  quadrature  to  calculate  the  far-held 
Green’s  functions.  The  near-helds  are  directly  calculated  using  an  adaptive  Clenshaw-Curtis 
quadrature  without  pole-extraction.  Results  of  numerical  simulations  are  presented  in  the  case 
of  a  graded-index  waveguide  and  an  atmospheric  gradient-layer  above  a  realistic  lossy  ground. 

1.  INTRODUCTION 

Electromagnetic  (EM)  modeling  of  refractive  and  terrain  effects  over  long  links  has  been  mod¬ 
eled  by  ray  tracing  [1,17,20]  or  by  using  the  parabolic/paraxial  wave  equation  [3,4,7].  These 
methods  are  approximations  to  the  wave  propagation  physics  of  Maxwell’s  equations,  which  can 
limit  their  applicability.  Other  approaches  have  included  the  moving- window  finite-difference  time- 
domain  (MWFDTD)  method  [14,16,21,22],  While  capturing  the  relevant  physics,  MWFDTD 
leaves  something  to  be  desired  with  respect  to  intuition  about  the  fields;  a  list  of  field  strength 
values  is  not  readily  understood  in  terms  of  propagation  mechanisms  such  as  direct  waves,  reflected 
waves,  interface/ground  waves,  and  guided  propagation  modes. 

On  the  other  hand,  the  theory  of  Sommerfeld  integrals  (Sis)  and  spectral  domain  Green’s  func¬ 
tions  (SDGFs)  in  multilayered  media  has  been  used  in  the  simulation  of  printed  circuit  board  (PCB) 
structures  such  as  transmission  lines  and  antennas  [8].  The  SDGF/SI  technique  expresses  the  EM 
fields  as  Sis,  which  must  be  evaluated  numerically  using  a  variety  of  methods  [2, 13, 15, 18].  With 
certain  approaches,  the  solutions  decompose  into  terms  that  directly  correspond  to  the  aforemen¬ 
tioned  propagation  mechanisms.  The  key  contribution  of  the  present  work  is  the  application  of 
these  rigorous  full- wave  techniques  to  long  range  radio  frequency  (RF)  propagation.  Quadrature 
methods  used  for  PCBs  are  not  directly  applicable  because  of  the  differences  in  length  scales  in¬ 
volved.  We  therefore  use  novel  asymptotic  quadrature  to  calculate  far-fields,  an  approach  that  has 
not,  to  the  best  of  our  knowledge,  been  attempted  for  RF  propagation  before. 

2.  ATMOSPHERIC  AND  EM  MODEL 

The  problem  geometry  is  a  non-magnetic,  multilayered  dielectric  material  bounded  by  two  half¬ 
spaces  with  interfaces  parallel  to  the  xy-plane  that  contains  a  z-directed,  time  harmonic  ( e~Jujt 
time  dependence)  dipole  radiator  at  a  height  z'  above  the  origin,  as  in  Figure  1.  Each  layer  is 
homogeneous  with  index  of  refraction  ni  =  where  sr(  is  the  relative  permittivity  of  the  £th 

medium  and  Eq  ~  8.85418782  x  1CU12  F/m  is  the  permittivity  of  free  space.  The  media  wavenumbers 
are  ke  =  n^ko,  ko  =  2irf  /co,  /  is  the  frequency  of  the  wave,  and  cq  =  299,  792,  458m/s  is  the  speed 
of  light  in  vacuum.  The  EM  fields  in  this  geometry  can  be  expressed  in  terms  of  only  the  z 
component  of  the  magnetic  vector  potential,  Az(p,  z,  z'),  where  p  is  the  xy-distance  away  from  the 
dipole  and  z  is  the  height  above  the  origin.  Az  obeys  a  forced  wave-equation  in  the  source  layer, 
and  an  unforced  one  in  the  others.  The  equation  reduces  to  the  3D  Helmholtz  equation  under 
harmonic  time  dependence,  which  is  indicated  in  Figure  1. 

The  spatial  domain  potential  can  be  shown  to  admit  a  spectral  (Sommerfeld)  integral  represen¬ 
tation  given  by 

POO 

Az(p,z,z')=  /  Az  (kp,z,z')  J0(kpp)kp  dkp,  (1) 

Jo 


Report  Documentation  Page 


Form  Approved 
OMB  No.  0704-0188 


Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 
VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 


1.  REPORT  DATE 

AUG  2013 


2.  REPORT  TYPE 

N/A 


3.  DATES  COVERED 


4.  TITLE  AND  SUBTITLE  5a.  CONTRACT  NUMBER 

A  Direct  Spectral  Domain  Method  for  Near-ground  Microwave  Radiation 

by  a  Vertical  Dipole  above  Earth  in  the  Presence  of  Atmospheric  5b.  grant  number 

Refractivity  5c.  program  element  number 

6.  AUTHOR(S)  5d.  PROIECT  NUMBER 

5e.  TASK  NUMBER 
5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES)  8.  PERFORMING  ORGANIZATION 

Georgia  Institute  of  Technology,  USA  report  number 

9.  SPONSORING/MONITORING  AGENCY  NAME(S )  AND  ADDRESS(ES )  10.  SPONSOR/MONITOR' S  ACRONYM(S) 

1 1 .  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release,  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

See  also  ADA586951.  Progress  in  Electromagnetics  Research  Symposium  (PIERS2013)  Held  in  Stockholm, 
Sweden  on  August  12-15,  2013.  AOARD  -  CSP-131024 

14.  ABSTRACT 

A  spectral  domain  method  is  presented  to  calculate  the  potential  of  a  vertical  dipole  in  a  multilayered 
medium  as  a  model  of  long-range  radio  propagation.  The  spectral  domain  Green’s  function  (SDGF)  for 
structures  with  up  to  hundreds  of  thousands  of  layers  is  calculated  using  an  e+cient  matrix  formulation, 
thus  enabling  simulation  of  continuously  strati  ed  media.  The  SDGF  is  sampled  to  perform 
pole/residue-extraction  using  contour  quadratures  and  to  perform  a  novel  asymptotic 
Filon-Clenshaw-Curtis  (FCC)  quadrature  to  calculate  the  far-  eld  Green’s  functions.  The  near-“elds  are 
directly  calculated  using  an  adaptive  Clenshaw-Curtis  quadrature  without  pole-extraction.  Results  of 
numerical  simulations  are  presented  in  the  case  of  a  graded-index  waveguide  and  an  atmospheric 
gradient-layer  above  a  realistic  lossy  ground. 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Georgia  Institute  of  Technology,  USA 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 


15.  SUBIECT  TERMS 


16.  SECURITY  CLASSIFICATION  OF: 


a.  REPORT 

unclassified 


b.  ABSTRACT 

unclassified 


c.  THIS  PAGE 

unclassified 


17.  LIMITATION  OF 

18.  NUMBER 

ABSTRACT 

OF  PAGES 

SAR 

5 

19a.  NAME  OF 
RESPONSIBLE  PERSON 


Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


Progress  In  Electromagnetics  Research  Symposium  Proceedings,  Stockholm,  Sweden,  Aug.  12-15,  2013  321 


(V"  +  k2L+l)  Az(L+1 )  —  0 

riL+ 1 

• 

• 

• 

• 

• 

• 

(V2  +  fc2)  Az,  =  -  r') 

ns 

(V2  +  k2)  A,3  =  0 

(V2  +  k l)  4*2  =  0 

n2 

(V2  +  k\)  42l  =  0 

n  i 

Figure  1:  A  dipole  radiates  inside  layers  of  dielectric  material.  There  are  L  layer  interfaces,  and  L  +  l  layers 
including  the  top  and  bottom  semi-infinite  half-spaces. 


where  Jo  is  the  Bessel  function  of  the  first  kind  and  order  zero  and  kp  is  the  radial  wavenumber  of 
a  spectral  component.  The  SDGF  A~  is  split  in  a  piecewise  manner  into  functions  in  each  layer, 
noted  as  Azp  for  £  =  1, 2, . . . ,  L  +  1.  These  functions  obey  a  ID  Helmholtz  equation  in  each  layer, 


d2A,u 

dz2 


+  kz£Az£ 


-HoIzO  S(Z4J  } ,  in  source  layer 
0,  in  other  layers, 


(2) 


where  Iz o  is  the  complex  magnitude  of  the  current  flowing  in  the  dipole,  /to  =  47r  x  10  '  H/m  is  the 
permeability  of  free  space,  and  in  each  layer  the  z  wavenumber  is  given  by  the  auxiliary  relationship 

kzt  =  \J k2  —  k2.  Equation  (2)  is  to  be  solved  in  all  layers  simultaneously  for  a  given  kp.  Then 
Equation  (1)  is  calculated  by  numerical  quadrature  over  various  values  of  kp. 

3.  SPECTRAL  DOMAIN  GREEN’S  FUNCTION  SOLVER 

The  general  solution  to  Equation  (2)  is 


7  -Poho 

Azt  =  / - ~ —  x 

**  8v r2 


+  R^e?ku(z  +  R£  e  ik‘Az  z>)  j  in  source  layer 


R+ejkzl(z~ze-i)  +  R~e-^z-^\ 


in  other  layers, 


(3) 


where  R±  are  generalized  reflection  coefficients  in  each  layer  that  have  to  be  fixed  by  boundary 
conditions.  The  boundary  conditions  are  that  the  tangential  EM  fields  have  to  be  continuous  at  the 
layer  interfaces,  and  that  there  are  no  incoming  waves  (Rf  =  =  0,  the  Sommerfeld  radiation 

condition).  Applying  these  to  Equation  (3)  reduces  the  problem  to  a  linear  system  of  equations  for 
the  Rf  coefficients,  the  exact  details  of  which  are  to  follow  in  a  future  publication.  The  system  of 
equations  is  sparse:  each  equation  involves  only  four  of  the  2 L  unknowns.  The  sparse  system  can 
be  assembled  into  a  2  Lx  2  L  pentadiagonal  matrix  with  at  most  (10L  —  6)  nonzero  entries  out  of  the 
full  4L2.  The  memory  requirements  for  large  numbers  of  layers  in  the  structure  are  not  prohibitive. 
Optimized  algorithms  exist  for  solving  pentadiagonal  systems;  this  leads  to  0(L )  time-complexity 
for  solving  for  Rrjr  ■  At  the  time  of  writing,  a  typical  computer  can  solve  systems  of  L  ~  100,000 
layers  in  ~  1  second.  Once  the  matrix  equation  is  solved  for  the  Rf,  the  SDGF  can  be  evaluated 
at  arbitrary  heights  z  using  Equation  (3).  The  matrix  is  formed  and  solved,  and  the  SDGF  is 
evaluated  at  multiple  heights  each  time  the  quadrature  routines  of  the  next  section  require  samples 
from  the  SDGF,  which  happens  during  both  pole-extraction  and  asymptotic  quadrature. 


4.  QUADRATURE  OF  THE  SPECTRAL  DOMAIN  GREEN’S  FUNCTIONS 

Sis  are  considered  to  be  difficult  to  numerically  integrate  due  to  their  slow  decay  and  oscilla¬ 
tions  [2, 12, 13, 15,  24],  The  slow  decay  is  from  poles  on  or  near  the  real  kp- axis  in  the  SDGF,  which 
degrade  the  ability  of  numerical  integration  algorithms  to  converge.  The  oscillations  are  from  the 
the  Bessel  function  Jo(kpp),  which  oscillates  more  rapidly  as  p  increases.  Conventional  quadratures 
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have  to  calculate  more  points  from  the  integrand/SDGF  as  the  range  increases  to  resolve  the  details 
of  the  Bessel  function.  Far-held  calculations  are  therefore  time-consuming  and  intensive.  In  this 
work,  the  slow  decay  and  rapid  oscillation  problems  are  handled  using  two  separate  techniques,  one 
well  established,  and  the  other  novel. 

4.1.  Handling  Slow-decay  with  Pole  Extraction 

Poles  and  residues  are  extracted  from  the  SDGF  and  analytically  integrated  using  the  residue 
theorem.  An  effective  search  algorithm  based  on  complex  contour  integrals,  which  automatically 
finds  the  pole  wavenumbers  and  residues,  has  been  described  in  the  literature  [9, 13].  In  the  present 
work,  the  poles  are  located  by  refining  user-provided  initial  guesses  with  a  direct  search  optimization 
and  then  contour  integrals  are  used  to  calculate  residues.  These  are  approximated  by  an  adaptive 
Gauss-Kronrod  quadrature  adapted  from  [19,  23].  The  residues  at  all  desired  heights  are  calculated 
in  parallel  by  the  chosen  Gauss-Kronrod  algorithm.  The  poles  are  known  to  come  in  opposing  pairs, 
±fcp,  so  the  residue  theorem  is  applied  to  Sis  of  the  form  /0°°  ,  2/,  2  Jo(kpp)kp  dkp  =  ^-H^\kvp), 

where  the  pole  is  at  kp  =  kp,  and  is  the  Hankel  function  of  the  first  kind  and  order  zero. 

4.2.  Handling  Oscillatory  Integrands  with  Filon-like  Asymptotic  Quadratures 

The  standard  approach  to  quadrature  of  oscillatory  integrands  is  to  sample  them  on  enough  points 
to  resolve  the  oscillations  and  then  approximate  the  integral  as  a  weighted  sum  of  the  sampled 
values.  For  the  present  case  this  would  mean  that  the  number  of  SDGF  samples  would  increase 
in  proportion  to  the  range  p.  Calculation  of  far-fields  would  require  more  computational  effort 
than  the  near-fields.  However,  there  are  several  families  of  asymptotic  quadratures  that  increase 
in  accuracy  as  the  oscillation  increases  [10, 11].  The  method  used  presently  is  the  Filon-Clenshaw- 
Curtis  (FCC)  rule  [5,  6].  Although  designed  for  complex  exponential  oscillation,  far-field  Sis  can  be 
adapted  for  the  FCC  in  the  following  way.  Since  we’re  interested  in  asymptotic  results,  the  Bessel 
function  is  replaced  by  its  asymptotic  expansion  in  Equation  (1),  giving 

Az(p,  z,  z')  «  \j-^-p  J  Az  (kp,  z,  z')  [e3kpP~^  +  e~3kpp+3^j  v /%  dkp,  (4) 

which  can  be  broken  up  into  two  integrals,  each  of  which  has  a  form  that  is  directly  amenable 
to  FCC  quadrature  after  the  semi-infinite  interval  is  truncated  to  a  large  finite  value.  Finite 
truncation  is  justified  physically  because  the  pole-extracted  SDGF  decays  exponentially  beyond 
the  largest  material  wavenumber  and  the  rapid  spatial  oscillations  of  large  wavenumbers  are  known 
to  contribute  only  to  the  near-field  singularity;  the  corresponding  waves  do  not  propagate  into  the 
far-held. 

5.  NUMERICAL  EXPERIMENTS 

The  method  presented  is  general,  and  could  be  applied  to  any  multilayered  problem.  As  a  proof- 
of-concept,  consider  a  PCB-type  example  in  which  a  copper  substrate  has  a  graded-index  material 
coated  on  it,  with  a  linear  gradation  of  10%  over  ten  wavelengths.  A  plot  of  the  extracted  guided 
mode  potentials  as  a  function  of  height  and  range  appears  in  the  top  panel  of  Figure  2.  There  are 
seven  poles  corresponding  to  guided  propagation  in  the  structure.  It  is  clear  from  the  plots  that 
the  structure  acts  as  a  waveguide,  trapping  the  energy  in  the  refractive  gradient.  In  this  example, 
the  source  dipole  is  3.1123  wavelengths  above  the  substrate.  The  lower  panel  of  Figure  2  is  a  plot 
of  ray  paths  that  propagate  in  the  same  gradient  over  a  reflective  ground  plane.  The  two  plots 
compare  favorably  in  terms  of  the  distance  scale  of  the  modal  “skip  zones.” 

Another  numerical  example  appears  in  Figure  3,  with  the  top  panel  representing  a  dipole  radi¬ 
ating  above  earth  with  er  =  15  and  a  conductivity  of  a  =  12  x  10-3  S/m  in  a  strong  gradient,  with 
n(z )  =  1  +  e~z/10  above  ground.  While  this  is  an  unusually  strong  gradient  for  atmospheric  cases, 
it  is  an  exaggerated  example  to  illustrate  the  effects  of  propagation  through  refractive  gradients  at 
ranges  that  can  be  easily  visualized.  For  more  realistic  gradients,  the  effects  are  apparent  farther 
away  from  the  dipole.  Because  the  FCC  quadrature  becomes  more  accurate  farther  away,  and  there 
is  no  additional  computational  cost  of  increasing  range,  there  is  no  inherent  problem  with  calculat¬ 
ing  extremely  far  fields,  only  one  of  visualizing  the  long  length  scales  in  the  kinds  of  plots  presented. 
Log-scaled  line  plots  can  be  appropriate,  but  we  find  the  visualization  of  the  figures  presented  more 
illustrative.  The  visualized  potential  is  calculated  by  the  sum  of  guided  modes  and  the  asymptotic 
FCC  quadrature  result  for  p  >  20 Aq  and  by  a  direct  Clenshaw-Curtis  quadrature  for  p  <  20 Aq.  The 
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Dipole  guided  modes  at  5.8  GHz  in  a  linear  gradient  above  copper 
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Figure  2:  A  dipole  radiates  in  a  layered  coating  on  a  copper  substrate  at  5.8  GHz.  The  potential  Az(p ,  z )  is 
visualized  on  a  decibel  scale. 
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Dipole  radiation  over  PEC  with  no  gradient 
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Figure  3:  A  dipole  radiates  in  two  similar  scenarios.  The  top  panel  is  using  the  SDGF /SI  method  in  an 
actual  scenario  and  the  bottom  panel  is  using  image  theory  in  an  idealized  scenario. 


fact  that  there  is  no  discernible  discontinuity  in  the  fields  at  p  =  20Ao  is  evidence  of  the  accuracy 
of  the  far-held  approximation  and  FCC  quadrature.  The  lower  panel  of  Figure  3  is  a  comparison 
with  a  simplified  environmental  model  that  can  be  solved  using  image-theory,  a  technique  that  is 
widely  used  by  communications  engineers  as  the  “two-ray”  model.  The  differences  between  the  two 
are  apparent  at  the  surface  level,  where  the  two-ray  model  can  be  said  to  break-down. 

6.  CONCLUSION 

A  direct  method  for  calculating  EM  fields  in  multilayered  media  based  on  SDGFs  and  quadrature 
of  Sis  has  been  presented.  Application  of  asymptotic  quadratures  to  SI  problems  is  novel  to  our 
knowledge.  Although  the  technique  is  conceived  of  as  a  method  for  simulating  long  range  radio 
propagation  near  ground  level,  it  is  applicable  to  general  EM  propagation  in  multilayered  media. 
Initial  testing  indicates  agreement  with  the  physics;  indeed,  the  calculated  fields  do  solve  Maxwell’s 
equations  and  boundary  conditions.  Future  directions  include  careful  validation  against  other 
models  and  measurements,  exploring  enhancements  to  speed  convergence,  formulating  the  solution 
for  horizontal  dipoles,  and  finally  handling  terrain  through  boundary-integral/scattering  methods. 
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